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Abstract. The formulation of a model for the evolution of the flow of a solid-liquid mixture (coal-water) 
in a horizontal pipeline with partial phase separation is the aim of this work. Problems of instabilities 
due to complex eigenvalues, observed in previous models, seem to be completely solved in the present 
model, in which we give the genesis of the different terms written in the equations, coming from the 
natural definition of mass and momentum balance, and the consequent proof of well-posedness of the 
obtained PDE system with boundary-Cauchy data. 

The model describes a three-layer flow. Most of the material is carried by the upper layer, while the 
bottom layer consists of an immobile sediment. The intermediate layer grows to a maximum thickness 
and has the role of regulating the mass exchange between the extreme layers. 

In the last section we present some simulations for a particular choice of flow regime, and boundary- 
Cauchy data, that were suggested by experimental results provided by Snamprogetti (Fano, Italy). 
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1. The physical problem 

In this paper we formulate a mathematical model for the flow in a horizontal pipeline of a solid-water 
mixture with a progressive phase separation due to gravity. 

The specific case we have in mind is the one of "dilute" coal-water suspensions and was suggested to 
us by Snamprogetti (Fano, Italy). The fresh mixture has a coal concentration of about 50% (in weight) 
and a particle size distribution centered at 0.205mm (< 1.25mm). The rheological properties of such a 
system are totally different from the ones of coal-water slurries that have been studied extensively in a 
number of papers (§, §, §, §, @, @, @). 

Indeed slurries, which can have coal concentrations up to 70% and are prepared using smaller particles, 
are stable at rest thanks to the action of chemical additives and exhibit partial sedimentation under 
stress, mainly because of the presence of impurities that are not affected by the additive and are no 
longer sustained by the slurry yield stress in dynamical conditions. 

The situation here is completely different, because the liquid and solid components tend to separate 
under gravity both at rest, and in dynamical conditions. 

(*) WORK PARTIALLY SUPPORTED BY THE ITALIAN CNR STRATEGIC PROJECT ON FLUID DYNAMICS AND BY 
THE ITALIAN MURST PROJECT ON MATHEMATICAL ANALYSIS OF PHASE TRANSITIONS 
(**) PRESENT ADDRESS 
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The approach we are going to follow is in the typical framework of stratified flows: in any transversal 
cross section of the pipe we have three layers, whose motion and composition are described by averaged 
quantities. 

The upper layer carries most of the material, while the bottom layer is immobile (a stationary deposit 
that should be avoided in practice). The middle layer has a density intermediate between the lighter, fast 
flowing, upper phase, and the sediment. We suppose that at the beginning just one phase exists, having 
the same properties of the material entering the pipeline. Therefore the system goes through a first stage 
in which there is no immobile layer, while the denser flowing layer grows up to a thickness A, depending 
on the pipe discharge. Once this value of the thickness is reached, the stationary sediment appears and 
the intermediate layer moves upwards, keeping its thickness A and its concentration, and transmitting 
solid material from the main phase to the sediment. 

The idea of the intermediate layer is inspired to a model that has succesfully explained sedimentation 
in slurries (see Q , [ fl6| ) , although as we said ,the general properties are quite different in the two cases. 

The scheme utilizing piecewise constant velocities over a cross section, although very convenient and 
frequently used in pipeline modeling, is already largely approximated, so that, taking the density constant 
in each layer, is not only a reasonable assumption, but it looks the only one really consistent with the 
general setting of the problem. In fact, the way density varies, over a cross section, can be measured 
by means of a Gamma-Densimeter, however the accuracy is not such that one could look for a more 
sophisticated model. 

Previous models in the literature describing flows of mixtures (see g, @, |]|, @, 0, Q, 
p2| , |24[ ] ) could not be adapted to our case for different reasons^ , so that the model we are going to 
illustrate is, as far as we know, an original contribution towards the study of non-steady suspensions 
flows with phase separation. 



2. The two-layer flow model 

Let us deal first with the initial stage of the flow in which we have just two layers. 
We denote by x the coordinate along the flow. At any section x = const., the two layers occupy two 
regions of areas Ai (the upper layer), A2 (the lower), separated by a straight line of length Si. The 



1 Problems of ill-posedness in []19| were pointed out (see Remark p.l| ) in pa] , while the stationary models present in 
the literature are not interesting in our case, since stationary flow is reached at such large dinstances from the entrance of 
the pipeline, that it is not observable in practice. In addition, it would correspond to the situation in which only water is 
transported by the main flow, which is of course of no practical use. 
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perimeters of the outer boundaries are Si, S2, respectively. In the stage we are considering the thickness 
of the lower layer is less than A (see fig. 2.1). 




\ 5; 

R 



h< A 




S 2 



Figure 2.1. Two- layer flow model, vertical section 



The flow model is based on mass and momentum conservation accounting for the exchange of liquid 
and solid between the two layers. 

The basic quantities entering the model are, besides A±, A%, 

a: solid volume fraction in layer 1; 
(3: solid volume fraction in layer 2; 
U: velocity of layer 1; 
V: velocity of layer 2. 

The solid volume transfer rate from layer 1 to layer 2 will be assumed to be a.A\ip, where ip is a positive 
constant. The corresponding quantity for the liquid component can be written (1 — a)A\Lp where (p turns 
out to be a function of a, (3, ip, as we shall see. 

2.1. Mass balance. Both the components are incompressible, so we have the following volume conser- 



vation laws: 



LAYER 1-solid 



(2.1) 



d t {aAx) + d x {aA x U) 



aA\ip 



LAYER 1-liquid 



(2.2) 



d t [{l-a)A 1 } + d x [(l-a)A 1 U] 



(1 - a)Anp 



LAYER 2-solid 



(2.3) 



d t ({3A 2 ) + d x {!3A 2 V) 



aA\ip 
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• LAYER 2-liquid 

(2.4) d t [(1 - f3)A 2 ] + 8 X [(1 - 0)A 2 V\ = (1 - a)A w 
Clearly A\ and A 2 are related by: 

(2.5) A X +A 2 =A 



Also the two velocities U,V are not independent, since summing up equations (2.1). (2.2), (2.3), (2.4), 
we obtain that the discharge Q along the pipe does not depend on x, as a natural consequence of the 
incompressibility of the mixture: 

(2.6) d x {A x U + A 2 V) = , 

so that Q may depend on t only, thus: 
(2.7) 



A 2 

The above equations are not enough to describe the evolution of the bottom layer. We assume that the 
solid concentration in it is constant, i.e.: 



(2.8) 



(3 = const. 



Thus we can divide eq. (^J) by and ( |2.4| ) by (1 — [3) and realise that 

a 1-/3 



(2.9) 



(p - 



1-a (i 



and that (2.3), (2.4) represent the same equation, as well as (|2.2|) is a consequence of (2.1), (2.3) and 



(2J 



Therefore we reformulate all the prevoius conservation laws using two equations only in the three 



unknowns a, A\, U. For instance we can replace A 2 and V in ( |2.3| ) using the expressions (2.5), (2.7) 
obtaining: 

a 



(2.10) d t A 1 + d x (A 1 U) 

and then we can modify ( |2.l|) as follows: 



Anl> 



(2.H) 



d t a + Ud x a = -a ( 1 - — ] ip 



We need one more equation which has to be obtained from the balance of momentum. 

Before deriving the missing equation, we define the average densities p\, p 2 in the two layers: 



(2.12) 



pi = ap s + (1 - ct)p n 



(2.13) 



Pi = fiPs + (1 - P)pv 
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(p s and p w are the densities of the solid and of the liquid components). 



While pi is constant, p\ varies. From (2.1) and (2.2), we obtain the global volume balance in A\\ 

(2.14) d t A 1 +d x (A 1 U) = ~A 1 i(> 

P 



after using (2.9). Similarly we have 
(2.15) 



d t A 2 + d x (A 2 V) = -Art , 



which is nothing but (2.5) divided by [3. 



Now, multiplying fl2.l| ) by p s , ( |2.2j ) by p w , and adding the results, we get: 



(2.16) 

and since the left hand side is 



d t {piAi) + d x (piAtU) = -A lP2 -ip , 



-pxAipip + Ai 



we deduce how pi varies along the flow: 
(2.17) D lPl 



dpi LT dp! 
dt dx 



i>{p2 - pi) 



where, here and in the following, D\ and denote the Lagrangian derivatives along the respective 
flows: 



1 dt^dx 







d 



D 2 = — + V— 



dt 



dx 



(2.18) 

As we are in the case a < (3 and p s > p w , and 

(2.19) p 2 - p x = 09 - a) ( Ps - Pw ) , 

we find out that (from ( [2.17 )) p\ is strictly decreasing along the flow, as long as a > 0. 

2.2. Momentum balance. Now we write the following pair of momentum balance equations per unit 
length: 

• LAYER 1 

(2.20) 

• LAYER 2 
(2.21) 



DMiPiU) = AiG - nSt T nSi + UDxipxA,) 



D 2 (A 2P2 V) = A 2 G - t 2 S 2 ± nSi - P2 UD 2 A X 



where G denotes the absolute value of the pressure gradient, t\, t 2 are the stresses at the respective wall 
portions, is the interfacial stress, with "+" , in ( p. 21 ) if U > V. We will come back to the form of such 
stresses later on . 
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Note that since DjA\ < for j = 1, 2 and D\p\ < (see ( [2.17 )), we have a momentum loss in (2.20) 
and a momentum gain in ( [2.21 ) . 



The genesis of the loss and gain terms is the following. Consider the motion of volume element A\dx 
of the upper layer. In the time clement dt it will transfer the mass element \Di(p\Ai)\dx dt to the other 
layer, with the corresponding (negative) momentum transfer UDi(piAi), per unit length and unit time. 

If we follow the motion of an element A2 dx of the lower layer, it will acquire the volume D2A2 dxdt 
in the interval time dt. The corresponding mass increase is obtained multiplying by P2 (constant) the 
volume accretion, and by U , in order to get the momentum variation. Remember that the matter coming 
from the upper layer, had a flow velocity U and that it is slowed down to V after passing the interface. 



We can rewrite ( 2.2C ) and ( 2.21 ) in a more expressive way: 

(2.22) p 1 AiD 1 U = AiG-nS 1 TnS i 

(2.23) p 2 A 2 D^V = A 2 G - t 2 S 2 ± nSi + p 2 (U - V)D 2 A 2 

Remark 2.1. The first equation says that the "donor" system obeys Newton's lau^ 

The second equation says that the accepting system experiences an extra inertia force due to the relative 

momentum of the transferred element. 

Note that, due to the relative motion of the layers, the mass element \D\(p\Ai)\dx dt is distributed 
among different elements of the lower layer. Similarly, the mass increase P2D2A2 dx dt is the sum of the 
contribution of different elements of the upper layer. Therefore (2.22), ( 2.23| ) should not be used to get 
the overall momentum balance. 



2.2.1. The stress terms. In ( 2.20| ), ( [2.2l| ) we take account of the wall and interface shear stress. Generally 



speaking, we introduce wall-specific shear stress (see jy): 
(2.24) r = ^pv 2 

where 

'Dv 



(2.25) A = c I 

' v 

is the friction factor (see jl2| , ffl) , v is the phase velocity, p is the average phase density, D is the equiv- 
alent hydraulic diameter, namely (W): 



2 Decompose the mass element piAi dx into (pi + D\p\ d£)(A\ + D\A\ dt)dx, which is still in layer 1 after dt and the 
complementary part —D\p\ dt(A\ + D\Ai dt)dx — p±DiAi dtdx (which is 0(dt)). The main portion moves according to 

(pi + Dipi dt)(Ai + D±Ai dt)DiU = external forces + interaction with the 0(dt) portion 



which gives ( 2.21 ) as dt — > 
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v a > Vh 


V a < V b 


V a ^ V b 




4A a 




-sfc— 


D b 


Sh 


Sh + Si 


Sh 



with usual notations, and the coefficients c, n can be evaluated, in the case of Newtonian fluids, for 
stratified flow, depending on the flow regime (laminar, turbulent) as: 



reg. 


c 


n 


Lam. 


16 


1 


Turb. 


0.046 


0.2 



In (2.25), v is the kinematic viscosity: 
(2.26) 



V 

v = — 



In the expression of the interface stress, v has to be interpreted as the relative velocity (U — V) = W of 
the two layers so that: 



(2.27) 

Here Ai and pi are the friction factor and the density of the faster phase. 



2.2.2. A model for the viscosity. The viscosity of the mixture is related to solid fraction. For this reason 
we select a particular behaviour of 77(a), considering, on one hand, the typical value of the viscosity 
of water (0.01gcm~ 1 s~ 1 ), corresponding to a = 0, and, on the other hand, the value of the viscosity 
measured for a suspension with solid fraction 0.5, typically 0.35gcm~ 1 s~ 1 (Snamprogetti). Since the 
bottom layer with solid fraction 0.7, is at rest, we impose that viscosity increases considerably beyond 
the value 0.5 for a. These considerations suggest the following choice of 77(a): 



(2.28) 



r/(a) = f] w [ ^- 

, '/to 



(see Fig. 2.2) where rj w is the viscosity of water and b is a parameter to be determined in such a way 
that — = 35, for a = 0.5, and 10 3 , for a = (3 (= 0.6). 

Remark 2.2. The selection above ofrj(a) is extrapolated from just two experimental data, so that it can 
be used to obtain just a qualitative, although quite realistic, description of the system. 
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1000 F 



100 r 



10 




Figure 2.2. Exponential model for the viscosity of the mixture 



3. The final equations 



Using the definition of D±, D% in (2.18), dividing (2.22) by A\ and (2.23) by A 2 and subtracting the 
second equation from the first, we obtain: 

(3.1) (A 2Pl + A lP2 ) d t U + (A 2Pl U + A lP2 V) d x U = P2 Q - ^nSx + t 2 S 2 t nSi + 1 



which is the third differential equation of the model (together with (2.11) (2.10)), in the three unknowns 



a, A\, U. Equation (2.7) gives 



V = V{a,A 1 ,U) 



while P is given. 



We can easily rewrite ( 2.11 ), ( 2.10 ), (3.1) in matrix form, introducing M, N, 3 x 3 matrices and O, 
F, two column vectors: 



(3.2) 

where 
(3.3) 

(3.4) 



ft 



a 
U 



V / 



/ = P2Q - j^nSi + r 2 S 2 t nSi f ^ + 1 



10 \ / u 

\[ =|01 \ N = U At 

A 2Pl +A lP2 I \00 A 2Pl U + A lP2 V 



so that: 
(3.5) 
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Md t fl + Nd x fl = F 
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Theorem 3.1 (Well-posedness) . The problem (| 
(3.6) 

is well-posed. 



with boundary- Cauchy data: 



n(o,t) = il° ; 
n(a;,0) = ft ; 



Proof. We begin by showing that (3.5) is hyperbolic (see ||, 0], (ITT)). M is invertible, so multiply (3.5) 
on both sides by M _1 and obtain the normal form of hyperbolic systems of PDE's 



(3.7) 

where now: 

(3.8) 

and 

(3.9) 

where 
(3.10) 



d t fl + Ld x fl = H 



L = M~ N = 



U 
u 







A 2Pl U + A lP2 V 
A 2 p\ + Aip 2 



H = M _1 F = 



h ) 



h = 



-, -— : — [piQ - ^r T iSi + t 2 s 2 t nSi (4^ 

A 2 pi + Aip 2 \ Ai \Ai 



As it can be easily seen, L has real positive eigenvalues: 
(3.11) 



Xi,X 2 = U 

. A 2Pl U + A lP2 V 

^3 — 



A 2 pi + Aip 2 

Moreover, as the three eigenvalues are positive and finite, the boundary-Cauchy data on the two axes 
(x = , t = 0), are given on time- like segments (see ||, [0), and the well-posedness of the problem is 
guaranted. □ 



Remark 3.1. Ill-posedness of previous models formulated for similar non-steady, multiphase flow, (see 
J15IJ had already been pointed out in different papers (see Jl^j , JTH , [p4|j where different methods to bypass 
ill-posedness were proposed, based on the introduction of surface tension terms. Such a procedure, which 
would be not applicable to our case, is not necessary in our model which is consistently formulated as a 
correct evolution problem in a natural way. 
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3.1. Non-dimensional variables. Before deriving the complete model with the introduction of the 
third layer, let us define a typical length-scale L c (e.g. the length of the pipeline, between two pumping 
stations) and a typical time-scale t so to use non-dimensional variables; the time-scale will be chosen to 
be jj^, with C/°the velocity of the fresh mixture at the entrance of the pipe. In particular L c will be the 
average distance between two pumping stations (~ lOOfcm), while U will be 2.0tos~ 1 in accord with the 
expected total flow rate (~ 1424m 3 /i _1 , (0.4m 3 s -1 )), with a section diameter of 0.5m. 



(3.12) 

Now we can use non-dimensional unknowns 



t x 



(3.13) 



Ol,2 



«1,2 = 



(U,V) 



A ' ~ x ~ U° 

We use also non-dimensional densities, scaled with the density of water p u 

Pl,2 



(3.14) 



Si, 



Pu 



non-dimensional cross-sections perimeters, scaled by the perimeter of the pipeline: 

S\2.i 



(3.15) 

non-dimensional specific shear stresses 
(3.16) 



°"l,2,: 



Ml, 2, 



Tl.2, 



where Ai is the friction factor of the wall shear-stress for the upper layer, evaluated at the boundary. 
We can rewrite the system (|3.7|) in a completely non dimensional form: 



(3.17) 

where: 

(3.18) 

with 

(3.19) 

and 



<9 r U + Ad^XJ = S 



U 




--ai*V 
P 



1 



a 2 5i + aid 2 



Ai 



o 2ttRL c 
A 



a-2 



P-l G\ + P202 T Pi hi (T 



ai 



«2 



«1 







( 







(3.20) 


A = 













\ 









ai 

SiajVi + S 2 aiV2 
Sia 2 + 8 2 ai 
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Remark 3.2. Here and in the following we considered constant the total volume discharge along the pipe, 
i.e. 



(3.21) 



Q = AU° 



4. The three-layer flow model 



Let us now introduce the third layer, seen, as we said in sec. ^, as a stationary deposit with constant 
solid volume fraction 7, and let us write the volume conservation equations, introducing the new transfer 
rates ip, Tp, of the solid and liquid components from layer 2 to layer 3. Then, the volume conservation 




A 



S3 

Figure 4.1. Three-layer flow model, vertical section 

equations will be: 

• LAYER 1-solid: 

(4.1) d t (aAi) + d x {aA 1 U) = ~aA^ 

• LAYER 1-liquid: 

(4.2) d t [(1 - a)4i] + d x [(1 - a)A x U] = -(1 - a)Ai<p 

• LAYER 2-solid 

(4.3) d t (pA 2 ) + d x (PA 2 V) = aA 1 (t/j - i>) 

• LAYER 2-liquid 

(4.4) d t [(1 - 0)A 2 ] + d x [(1 - 0)A 2 V] = (1 - a)A x (if - Jp) 
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• LAYER 3-solid 
(4.5) 

• LAYER 3-liquid 
(4.6) 



d t (jA 3 ) = aAx-ijj 



d t [{l-j)A 3 ] = {l-a)A x Tp 



Remark 4.1. The reason why the quantity of matter passing from the middle layer to the stationary 



deposit is proportional to a.A\, and not to the corresponding factor of the middle layer PA 2 , in (4.5), is 
due to the fact that we suppose the thickness A of the middle layer constant . From this it follows that 
the transfer rate across the intermediate layer is driven by the trasfer of matter between the upper and 
the middle layer. 

In this way, the mass transfer to the stationary deposit stops, as well as the one from the upper to 
middle layer, only when the upper layer is completely empty of solid (a = 0). This is the condition of 
complete phase separation and consequent steady flow, which is reached asimptotically by our system, 
although in a region of not practical interest, since it is far from the operating conditions of a plant. 



Now we have 

(4.7) 

and 

(4.8) 



A = A t + A 2 + A 3 



/3, 7 = const. 



Summing up equations from (LI) to ( |4.6| ), we still have the total volume conservation: 
(4.9) d x (A 1 U + A 2 V) =0 

which means again: 
(4.10) 



A-2 



Dividing (4.5) by 7 and (4.6) by (1 — 7), after using (4.S), and equating the results, we get: 
(4.11) Tp 



a 1 - 7- 
-ip 



1 — a 7 



which replaces (2.9) 



Dividing (L3) by (3 and (L4) by (1 — (3), and equating the results, we have now: 
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Note that, being the thickness of the middle layer constant, we can express 
(4.13) A 2 = A 2 (A 3l A) 

from which it follows: 

8A 2 dA 2 dA 3 



(4.14) 



d(x, t) dA 3 d{x, t) 



Remark 4.2 (On the thickness of the middle and bottom layer). To simplify the equation and without 
losing physical meaning, we will suppose that the thickness of the sedimentation layers can be considered 
small enough to have (see Fig. \4-j ): 

{pi 3 

(4.15) sinifii ~ ifi — , i = 1,2 



(> 



With this hypothesis, we will have: 



(4.16) A 2 ~\r 2 { v \- v I) 



3 



and 



(4.17) A 3 ~ \R 2 V \ 



3 



while tpi and ip 2 will be related by: 



(4.18) A ~ ±R(cpl - V i) 

This assumption can be made without any loss of physical meaning since we know that the presence of 
the stationary deposit is a damage for the correct use of the plant. Therefore, in practical applications in 
plant, its thickness must be kept below a small fraction of the pipe radius. 



From (4.7) we get 

dA 3 dA 1 dA 2 



^■ 19 ^ d(x,t) d(x,t) d(x,t) 

Using (4.14) in ( 4.19Q , we have then 

dA 3 1 dAt 



(4.20) 



Now substitute (00) in (|4.5|): 



d{x,t) dA 2 d(x,t) 

dA 3 



(4.21) ^=~ 2 ^4-¥ i 

a -. oa 2 A\ 

dAl 



which completes the relation between if and ip in (4.12) 



14 



ALESSANDRO SPERANZA 



Remark 4.3. We will actually suppose that the transfer of matter from the middle layer to stationary 
deposit starts only when the thickness of the middle layer has reached the constant value A. We will 
suppose in particular that: 

7 1 dtAi. 



a dA 2 Ay 
OA, 



-T(h) 



where T(h) is the step function of the thickness h of the middle layer: 



Using (4.8) we can rewrite (4.3) as: 
(4.23) 



d t A 2 + d x {A 2 V) = -,Ay ty-ty 



and using ( ^14J ) together with ( t4.20| ), ( |4.2l| ) and the definition of V, ( f4-10|) , we get: 

8A 2 , 7 

(4.24) 



9M dA 2 dtAl + d ^ A ^ = 



1 + 



dA 3 



Now we divide (4.1) by a and subtract (4.24) from the resulting equation: 
Ay 



(4.25) 

that we can rewrite 
(4.26) 



-d t a + 



dAz [3 



1 - 



V 



1 



dA 2 
dA 3 ) 



d t Ay + ^d x a = -Ay(l-^)r 



1_1 

dta ~ a ~ 8A 2 ^a! + U9xa = ~ a i 1 ~ ij ' ' 



1 



dA 3 



We will take (4.26) and (4.24) as the the first two equations of the final model in the unknowns a, Ay, U. 



Remark 4.4. Equation (4.26) can be rewritten, using the usual definition of the derivatives along the 
flow, 



(4.27) 



1-1 

Dya = -a[l--U + a--g^ — 



1 



dA 3 



and, since dtAy < 0, a < (3 and (3 < a is strictly decreasing along the flow. 

4.1. Momentum balance. In the following we will use the same notations as in sec. [2] with the changes 



due to the introduction of the stationary deposit, as described in Fig. 4.1 



The momentum exchange will be now between the upper and the middle layer, with exactly the same 
terms as before and a term of momentum loss of the middle layer, in favour of the stationary deposit. 
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The momentum transferred to the stationary deposit is completely absorbed by the pipe wall. With this 
in mind, let us write down the two equations of momentum balance for the layers in motion. 

• LAYER 1: 



(4.28) 

• LAYER 2: 



Di(AipiU) = A±G — t\S\ =f t 12 S 12 + UDi(ftiAi) 



(4.29) 



D 2 (A 2P2 V) — A 2 G — t 2 S 2 ± T12S12 — T23S23 — p 2 UD 2 A\ — p 2 VD 2 A 3 



Using Q4.2CD , ( |4.10| ) in ( |4.29[ ), dividing Q4.28Q by and ( f4.29| ) by A 2 , and subtracting the first from 
the second, we get the third equation of the complete model: 



(4.30) 



P \ A d t A x + (A 2Pl + A lP2 ) d t U + 92 d x A x + (A 2Pl U + A lP2 V) d x U 
dA 3 + dA 3 

= ~~^~ Tx ^ 1 + T2 ^ 2 ^ T 12^12 + 1 J + T 23S 23 



Before going further with the calculations, let us reduce equations ( 4.26 ), ( 4.24| ), ( 4.30 ) to non- 
dimensional form, with the same scale factors as before: 



(4.31) 



ct ( a\ n 

f- o T ai + videa — —all — —) rip 

} 02 +1 ai \ 0) 



da 
da 3 



(4.32) 



(4.33) 



da 2 7 

g^-o T ai + (aivi) = —ai—t"ip 



da 3 



l-l ^ 
^ — 2V2 d T ai + (a 2 Si + ai5 2 ) d T vi + 

1+ da 3 

1-1 

^2 2 a , 1 x x \ » 
g^- -J a^ai + (a 2 0iwi + 0102^2) 0£Vi = 



1 + 

2L C ( a 2 ( a 2 

Ai — — fixat + P2 a 2 =F hi P \ 2 o\ 2 + P23 o 23 

ti \ a\ \a\ ' 



The boundary-Cauchy data, will be a constant vector: 



(4.34) 



u" = ! 1 
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Theorem 4.1 (Well-posedness) . The system of PDE (4.31), (4.32), (4.33) with boundary- Cauchy data: 

f4 35) fU(0,r)=U o ; 
(4 ' 35) i U(£,0) = U°; 



is well-posed in the range of validity of the assumption (4.15) 



Proof. The system of equations (4.31) ( 4.32| ) (4.33) can be reduced again to normal matrix form 



(4.36) 

where, as usual: 



<9 T U + Ld s U = S 



(4.37) 



U 



Q 

ai 

Vi 



1 ^1 _i_ 1 P 
V da 3 (3 J 
oa 2 



t°4> 



da 3 



1 



da 2 7 P 
da 3 (3 
s 



ai-iV 



with 
(4.38) 



1 a ai5 2 v 2 ,o , 
t tp- 



da 2 7 (3 a 2 5i + aifo 
<9a 3 /3 

Ax 2X C / a 2 



a 2 5i + ai5 2 R \ ai 
and where the matrix L is now: 

/ 



a 2 



— Ml "! + M202 T I — + 1 ) M12C-12 + ^230-23 



7 



(4.39) 



1 



/3 av 1 (3 



1 \ 





V 



da 2 7 ai <9a 2 7 

da 3 /3 9a 3 [3 

1 + P 1+ |02 

gag L gag 

7 , £^2 1 7 , £^2 

/3 9a3 (3 da 3 



a 1 



32 



33 



with: 
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It is easy to see that, defined a function 



(4.41) 



2 





da 2 7 ' 
da 3 (3 



the matrix L can be rewritten 

(4.42) L = A + eC 

where A is the matrix of the two-layer model 



(4.43) 



and C is 



(4.44) 



.4 





a i 



\ 



h02Vl + 5 2 CLiV 2 

8\a 2 + 82CL1 ' 



( 





c 



ai 

-Vi 



a \ 

-ai 



-k 



V 



vi - v T 



dd2 + Z\ 

da 3 13 



\ 



1 



—a\k 



J 



da 3 } 

so (see Remark |4.6| ), for small e, the introduction of the stationary deposit can be seen as a small 
perturbation of the two-layer flow model. 



We want to show that the equation ( 4.36 ) is hyperbolic. The first eigenvalue of L is easily seen to be 
v\, the other two are the zeros of the polynomial: 

(4.45) p L (A) - (Z 22 - A) (Z 33 - A) - l 32 l 23 

To show that the eigenvalues are real, we show that: 



(4.46) 



(^22 + hzf — 4(^22^33 — hzh 2 ) > 



in the hypothesis of small stationary deposit. 
This can be rewritten: 



(4.47) 

which means: 

(4.48) A L 2 
where we used: 
(4.49) 

the third eigenvalue of the matrix A. 



Az, 2 = (J22 - hs) 2 + 4^32 > 



-i 2 



(1 — e)v\ — A3 + ea\k — 4ea±k {v\ (1 — e) — v 2 ) 



A _ a 1 8 2 V2 + a 2 8ivi 
aiS 2 + a 2 5i 
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Let: 
(4.50) 
and 
(4.51) 
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ai<5 2 



a\82 + 0,261 



7-/3 



Now we can rewrite eq. (4.48): 
(4.52) 
Calling 



(i-r e ) 2 + 4^ 

UJ 



v\ - 2ujv x [(1 - re) (uj - e) + 2eT(l - e)] v 2 + vf(uj - ef 



(4.53) e = Te, 

( 4.52j ) can be rewritten in the more compact way: 



(4.54) A L ' = vf(uj-e) 

where 

(4.55) 



(1-e) 2 +4- 



u l - 2 



uj — e 



(1-3 + 



2e(l - e) 



— e 



tt + 1 



U2 
Ul 



Therefore, the zeros of A 2 correspond to the zeros of a polynomial in u of the form: 



that has solutions for 



"1,2 



au A - 2 



few + 1 



e\ b±VV TZr a 



uj J a 

After some algebra we will see that these solutions are in fact out of the region in wich vi varies. 
Neglecting the e 3 terms (see Rem. |4.6| ), we find that: 



(4.56) 



CI ( i+i 



So, neglecting the e 2 terms in & and a, we end up with: 



(4.57) 



"1,2 



uj uj V VI 



1 + 2e I 1 



where e, e << 1 (see Rem. 4.6). 

From the definition of uj, we see that: 

(4.58) 



1 -uj 



a2<5i + a\52 
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which is of order ipf — Lp\ (see eq. (4.16)), so that in the expression (4.57), the difference between the 



two zeros u\ and is small in comparison with u\, u 2l which are both close to 

a 2 Si 



(4.59) 



1-e 1 



the difference being of the order of <^2fi 



2/3 



From the last we can actually appreciate that u has to be close to 1, to make negative, conditions 



that will never be reached in the hypothesis ( 4.15 ) (see Remark 4.5 ) 



□ 



Remark 4.5 (On the validity of Theorem 4.1). It is clear that the interval in which A^(u) is negative, 



is very close to u — 1, conditions in which we are not interested (see Rem. l^.l ), corresponding to large 



phase separation (see Rem. 4-0> negligible in practical applications. 

From asymptotic study (see App. we know in fact that v%, after a jump-like behaviour for f « 1, 
starts increasing, pushed by the faster, less dense, upper phase, but remaining considerably smaller than 
V\, also increasing. 



Remark 4.6 (On e). It is easy to see that, for (23 << 1, 

¥2 



(4.60) 



<p 2 + r[<P2 + 2 



1/2 



with tf2, the angle subtended by 03 (see Fig. 4-1). 
Therefore e << 1 if 

r2 



(4.61) 



r 2 



R 



172 « 1 > 



which defines what we mean by a thin stationary deposit. 



Remark 4.7. The case in which 03 is not small, like, e.g., in the case of large phase separation, should 
be treated differently and will not be considered in this model. For example in problems of restart, the 
pressure gradient, that we considered the same for the two layers in motion, is different at different levels, 
during a transient phase. 

5. ON THE BOUNDARY- CAUCHY DATA 



In sections H and 4.1 we introduced the boundary-Cauchy data. In particular, in theorems 3.1 and 



4.1 we imposed constant data on both axes (£ = 0, r = 0). 
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It has to be noted that this choice was not the only possible. We chose that at the entrance the 
condition of the mixture is, at every time f, the same as it was in the whole pipe at time t = 0. 

Note that the regularity of the boundary-Cauchy data influences the regularity of solution (see JllJ , 




Theorem 2.1, pag 71). In our case, being the functions L and S in (4.36) regular and choosing constant 
the boundary-Cauchy data, the solution will be regular too. 
Our choice will be therefore: 

(5.1) 

This corresponds to the choice we made of having at time t — and at any t, at the entering cross section, 
the pipe full of matter with solid fraction a and velocity U , namely: 

a = 0.5 A\ = A U° = 2m/ s 

5.1. On V°. Some remarks have to be made for the value of the velocity of the middle layer at the 
entrance. 



From eq. (2.7) we see that in fact V is defined only for A 2 ^ 0. However this is not true in our 
conditions, for x — 0. We can calculate the value V° as: 

(5.2) V° = lim Q ^-^ U 

x^O Ao 



To do this we have to take into account the different terms in the source function in eq. (2.22) and 



( 2.23 ), evaluated in the vicinity of x = 0. All the terms of the right-hand side vanish, linearly with A 2 
or with Si, as A 2 — > 0, apart from the wall shear-stress term T2S2 vanishing only in turbulent regime. In 
this case, V° turns out to be equal to U°, so the second layer appears with the same velocity as the fresh 
mixture. 

On the other hand, choosing the laminar regime, we get V° = 0, which is however consistent with the 
model only if ip is taken no longer constant, but a function of x, vanishig at x = 0. 

We put in the Appendix the study of the asymptotics near the origin, in which we distinguish the 
different cases of flow regime. 

In the numerical simulations we chose the turbulent-turbulent regime because it fits more naturally 
the model with tp constant. 

A possible genaralization could be the introduction of a transition from one regime to another in the 
first segment of the pipe. 
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6. The numerical results 
To evaluate the solution of eq. ( 4.36| ) we elaborated a numerical code, based on an explicit approxi- 



mation method . Eq. ( 4.36 ) is discretised as 
(6.1) 



, ri , Trir , (u\i] \xj] — u\i] \xj-i\) 

u[i + l][ Xj ] - u[{\[ Xj } = At \ L[i][xj] 3 Ax l 



H\i][ Xj 



where we indicate with u[z][xj] the solution of (4.36) evaluted at i — th time step, in the j — th cell of 
the spatial domain. 

This finite differences method of approximation explicit in time (Euler explicit) with a judicious choice 
of the dimension of the cells of the space-time grid, gives good confidence of stability and quite acceptable 
approximation, and allows the simulations to be run on a normal PC. In the numerical code we inserted 
a control for the dimension of the space-time grid, in order to satisfy the Courant-Friedrichs-Lewy^j 
condition (see UM) during the whole run. 



The solution of equation (4.36) can be plotted at different time steps as a function of the normalized 
distance from the entrance £ (= x/L c ). We plot, in the graphs below different quantities: a, V\ and 



1 — a±, solutions of ( 4.36 ), as well as some important ones, related to the solution vector (U), namely 
V2, a 3 . 



In Fig. 6.1 we show the result of a simulation of the evolution of the system in 100000 sec, as we 
said, in turbulent-turbulent regime. We see the formation of the stationary deposit after the achievement 
of the maximum thickness of the middle layer. Here we chose U° = 0.5m/s and L c = lOOfcm, with a 
maximum thickness A = 5cm. 

In similar conditions, but with a maximum thickness A = 8cm, we do not see any stationary deposit 



(see Fig. |6_J). 

With the initial velocity U° = 2m / s and A = 8cm, again, no stationary deposit is present on a lenght 



of lOOfcm (see Fig. 6.3), while the stationary deposit is visible in similar conditions, but with A = 5cm 



(See Fig. 6.4). From the above simulations we see the dependence of the deposit upon the thickness, 
A, of the middle layer and upon the input velocity U . The stationary deposit is, as we expect, thinner 



3 The condition states that, given the eigenvalues of the characteristic matrix, \ p , and being h (along x) and k (along 
t) the dimensions of the cells in the space-time grid, necessary condition for the convergence of an approximate method for 
PDE systems, is that: 

< 1 



In the case of the complete system (see sec. g), we calculate numerically the eigenvalues of the characteristic matrix and 
we reject the time steps at which the CFL condition is not satisfied. For a detailed discussion on the convergence of this 
approximation method, see 
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Figure 6.1. Numerical simulation of the three layers flow model for half a characteristic 
time (t° = 200000sec), A = 5cm 
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0.8 
0.6 
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0.2 




Figure 6.2. Numerical simulation of the three layers flow model for half a characteristic 
time (t° = 200000sec), A = 8cm 



in the case of a larger input velocity and disappears completely in figures 3.2 and 6.2 for large value of 
A. 

Note the jump-like behaviour of V2 near the origin which simulates the real behaviour V2 ~ x ' 2 (see 
App. g). 
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Figure 6.3. Numerical simulation of the three layers flow model for two characteristic 
times (t° = 50000sec), A = 8cm 
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Figure 6.4. Numerical simulation of the three layers flow model for two characteristic 
times (t° = 50000sec), A = 5cm 



We have also the convergence, observed for bigger values of the lenght-scale (out of the physical lenght 
of 100/cto) to an asymptotic solution of steady flow, corresponding to the condition of complete phase 
separation. Our model converges then to the steady models found in the literature, in B, pOf, p2[. 



24 alessandro speranza 

7. Acknowledgements 

This project was suggested by Snamprogctti (Fano, Italy), a leading company in the field of fuels 
pipelining. I am thankful to the researching staff of Snamprogetti, for providing most of the bibliographic 
material, and for the substantial contribution to the progress of this work. I thank also Dr. A.Mancini, 
from Universita di Milano, for his constant assistance in the development of the numerical simulations 
and Prof. A.Fasano from Universita di Fircnzc for his advise. 



SUSPENSION FLOWS IN PIPELINES WITH PARTIAL PHASE SEPARATION 25 

Appendix A. Asymptotic behaviour 

We try here to evaluate the behaviour of the solutions of the system of PDE' s we obtained, in the 
vicinity of the origin. 

Here we neglect the dependence of the solution upon t, situation that, as we observed in the develop- 
ment of the numerical simulations, can be actually reached in a time of the same order of the characteristic 
time t°. 

Let us suppose then that, for small x, the solution vector fl has the expansion: 



(A.l) 



Ax = A — kx m - k Q x A 2 = kx m + k x 

U = U° + h lX mi + hi°x V = V° + h 2 x m2 + h 2 °x 



with 771, mi, m 2 < 1. 

From eq. ( 2 . 1 d| ) , we get 



(A.2) 

From the conservation of total volume: 
(A.3) 



d x (AiU) =s -jAi; 



d x (AiU) = -d x {A 2 V) 



Using (A.l) in (A.2), equating the terms of the same power, we get: 



(A.4) 



, Ah^ - k U° = ~A%^ 
hi = k = 

Ahi° - k U° = -A—*P 
m = mx\ (3 

Ah x ~U°k = 



For m + mi > 1, while the case m + mi < 1 is not possible. 
If, in particular, 



(A.5) 



m + 777l = 1 



these conditions change a little in the case m — mi, while nothing changes in the other cases: 

n° 

, -khi + Ah^ - k U° = -A — ip 
(A.6) m = mi \ i^i u p y 

Ahi -U°k = 



Using now (A.3) in the same way: 



(A.7) 



k V° = A—il> 
m + m% 1 \ (3 

kV° = 



777 + 777 2 



I kh 2 + k Q V° = A?-i> 

\ fJ 



kV° = 
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Note that, in the case m + m 2 ^ 1, it cannot be V a = 0, while, with m + m 2 — 1 both, V° — and 
V° 7^ 0, are acceptable. 

We can do the same with the two equations of momentum balance fl2.22 ) ( 2.23 ); from the first we 
get: 

(A.8) Pl °AU° (m^!!""^ 1 + /i! ) ~ AG - Sn° - nS, 

Then it has to be: 



(A.9) 



hi = 



from which it follows (see (A.4) and ( |A.q )) 

(A.10) k = 

in all the cases. 

Moreover, as TiSi is small near the origin, 

(A.ll) //, 



o AG - St? 



Pl °AU° 

Then the behaviour of U around the origin depends on the pressure gradient. 



From ( p.23| ) we get: 

P2k x(V° + h° 2 x)(h° 2 + m 2 h 2 x m *- 1 ) ~ 

~ k xG - t 2 S 2 + r t S l + p 2 (U° -V°- h 2 x" 12 + [h\ - h° 2 )x) (V° + h 2 x m2 + h° 2 x)k 

(A.12) 

All the terms vanish in x — 0, apart from 

P2 (U° ~V°)V°k Q 

then, it has to be either 

y° = o 

or 

V° = U° 

In the case V° = U°, equating the linear terms we get also: 

< A - 13 > "5 = ^ + f 



and, as k is zero, from (A.4) and (A.6) we have 

(A.14) k = A ^0 + cP^ 
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Now let us call tp the angle subtended by A 2 , (see fig. 2.1). 
It can be seen that, for small <p, 

(A.15) r 2 S 2 ~ uj 2 V 2 - n ^- 2n 

where lu 2 is a positive constant, and n is the power appearing in the friction factor, (see sec. [2.2.1 ), 

depending on the flow regime. 

Note that, to prevent this term from becoming singular in the origin, we need V — > as x — > 0, if 
1 

U> 2- 



We can rewrite ( A.15 ) 
(A.16) t 2 S 2 ~ w 2 / V^- n ar i: T I 

where 

In the case n < — , it has to be then, 
- 2 



from which it follows that 
(A.17) 



1 - 2n 
m 2 = — 7: — 



LU 2 'U Ql ~ n 



p 2 k (m 2 + 1) 

So we have A 12 and U, linear in x, V ~ [/° + /i"^ + h 2 x^~ , with /12 negative. 
The case = is actually not consistent with the assumption ip — const. 

Let us see how we have to modify it in order to make this behaviour (suggested by the choice of laminar 
regime, n = 1) acceptable. 
The dominant terms are, 
• on the left-hand side: 



p 2 k h 2 h 2 x 



on the right-hand side: 



C)2 l+i 



T 2 S 2 ~ uj' 2 h 2 x 
nSi ~ WjJ7° x 
P2 U°k h 2 x m2 
p 2 k h 2 2 x 2m2 

2 

Then, if h 2 ^ 0, the t 2 S 2 term is balanced only if m 2 — — . This requires also: 
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(A.19) 

We said also: 

which means 
(A.20) 
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h° = u° 



d x (A 2 V) ~ A— V> 



koh 2 (m 2 + l)x" 12 + 2k h 2 °x ~ A— ^> 



The matching above requires a different choice for ip. For example we could say that, 



ip(x) = 



l/> X < Xq 

(A.21) 1 

In this way it has to be h 2 = and 
(A.22) 

and also V is linear around the origin; otherwise, if h 2 ^ 0, we must have: 

(A.23) V ^ V^™ 2 

and: 




(A.24) 



/l2 = 



a 

k Q (m 2 + 1) 



which is consistent with (A. 18), only if we take: 



(A.25) 



Wi k (m 2 + 1) 



oj 2 'U^ 



A- 



Now in (A.4) we would have 
(A.26) 

and a,A\, U, would stay linear. 



fe 



Ah\ 
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